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C/5 1 Abstract 

In this article, we propose a new numerical approach to high-dimensional par- 
tial differential equations (PDEs) arising in the valuation of exotic derivatives on 
the LIBOR curve. The proposed method is adapted from [20] and uses principal 
component analysis (PCA) of the underlying process in combination with a Taylor 
expansion of the value function into solutions to low-dimensional PDEs. The approx- 
imation is related to anchored analysis of variance (ANOVA) decompositions and is 
expected to be accurate whenever the covariance matrix has one or few dominating 
eigenvalues. A main purpose of the present article is to give a careful analysis of 
the numerical accuracy and computational complexity compared to state-of-the-art 
Monte Carlo methods on the example of Bermudan swaptions and Ratchet floors, 
which are considered difficult benchmark problems. We are able to demonstrate that 
for problems with medium to high dimensionality and moderate time horizons the 
presented PDE method delivers results comparable in accuracy to the MC methods 
considered here in similar or (often significantly) faster runtime. 

1 Introduction 

In most common models, the values of financial derivatives are equivalently characterised 
as the expected value of a payoff functional under some stochastic process or the solution 
of an associated partial (integro-) differential equation. The two dominant classes of nu- 
merical methods in derivative pricing are therefore Monte Carlo methods (see, eg, [B]) for 
estimating the expectation via simulation and discretisation methods (see, eg, [UI22]) for 
approximating the solution to the respective PDE (where we include lattice, spectral and 
Fourier methods in the latter group for the properties we shall discuss now). Simulation 
methods are well suited to track path-dependent quantities which determine the payoff of 
exotic derivatives, and scale favourably with the dimension of the process. However, the 
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convergence in the number of samples is slow and they require additional approximations 
to early exercise strategies. Conversely, conventional PDE discretisation methods incorpo- 
rate early exercise features easily and allow fast convergence in the number of nodes used 
in each direction, which makes them very efficient for low- dimensional problems, but they 
become intractable as the dimensionality increases. 

The effort to solve N- dimensional PDEs numerically with standard grid-based methods 
grows exponentially with N and even more sophisticated PDE methods tailored to high- 
dimensional approximation, such as those based on sparse grids, are typically not able to 
deal with practical problems where N exceeds about five to eight, see [9l [TUl [T3J, [20] . Given 
especially the advantages in dealing with early exercise, it would be not only of academic 
interest but also practically very relevant to be able to solve generic high- dimensional 
derivative pricing problems with PDE methods. 

In this paper, we adapt an approach from [2D] that computes an approximate solution 
of an iV-dimensional PDE by solving 0(N P ) PDEs of maximal dimension d <C N. In 
fact, we will see that p = 1 and d = 2 is usually sufficient for practically adequate accu- 
racy. The underlying principle of this and related approaches is an anchored ANOVA-type 
decomposition (see [IB]) of a solution u(z), z £ M. N , into 

N N 

u(z) = u (a) + ^Mj(a; z{) + ^ Ui,j{a; z u Zj) + . . . + Wi,...,jv(a; zi, ...,z N ) 

1=1 i,j = l 

i < j 

u v (a;z v ), (1) 

v C {1,...,JV} 

where we associate u% with Uq, with m etc. The terms on the right-hand side each 
only depend on a subset of the coordinates, z v = (z^, . . . , z^ v ,), and a chosen 'anchor' 
a = (ai, . . . , a at)- This has been successfully applied to quadrature problems from finance 
in [5] , and its relation to the PDE expansions in [20] , which form the basis for the present 
work, is highlighted in [19] and [2T] . 

Key to the efficiency of this approximation as a numerical method is that the size of u v 
decays rapidly with increasing \v\. This can be achieved by a coordinate transformation 
of the underlying stochastic process and of the corresponding forward or backward PDE. 
Optimal linear transformations taking into account the payoff function are analysed in 
[TTj . while here we consider the principal components of the covariance matrix E of the 
Brownian driver of the process. The accuracy of the approximate solution obtained by 
truncating ([Tj) after a small number of terms with small \v\ then depends on the (relative) 
sizes of the eigenvalues A« of S, 1 < i < N. This will be motivated in Section [2] by 
expanding the value function in Aj. We follow here [20], who first introduced this idea for 
vanilla basket options. 

In this article, we demonstrate the wider applicability in situations where no closed- 
form solution is known and accurate Monte Carlo estimates are difficult to obtain. A 
prime candidate for using this technique in practice is the LIBOR market model for the 
joint evolution of LIBOR rates with different tenors. To value path-dependent products 
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such as TARNs (Targeted Accrual Redemption Notes), Snowballs or Ratchets, and early 
exercise options such as Bermudan swaptions, indeed the whole yield curve has to be taken 
into consideration, which makes the problem genuinely high-dimensional for long enough 
maturities. The PCA-ANOVA-based PDE approach presented here is very well suited to 
this setting even in high dimensions, because LIBORs with similar tenors are strongly 
correlated, such that one observes a fast decay of the eigenvalues, as is seen from Fig. [TJ 
in Section |H On the example of Bermudan swaptions, even when including the first order 
terms with |u| < 1 alone, only a mild loss of accuracy is observed as the dimensionality, 
determined by the number of LIBORs considered, ranges up to 50-60. This deterioration 
appears to be an effect mostly of the time to maturity and for longer running contracts, 
the higher order terms in ([1]) become more relevant. Putting this in the context of the 
commonly used Monte Carlo method presented in [2J, the necessary restriction of the class 
of exercise strategies there produces a gap between lower and (dual) upper bounds which 
widens as the maturity increases. The accuracy of these Monte Carlo results is comparable 
with the PDE ones, which are obtained in a small fraction of the computational time. 
Additionally, the expansion (TjQ) implicitly defines a systematic accuracy improvement and 
is relatively straightforward to implement. 

This paper is organised as follows: Section |2J introduces the PCA approach, and Sec- 
tion |3] discusses its relation to anchored ANOVA decompositions. In Section |4] we apply 
the approach to the LIBOR Market Model, and show numerical results for two LIBOR 
derivatives, Bermudan swaptions and Ratchet floors, in Section [5j Section |6] summarizes 
the results and discusses extensions. 

2 A PCA-based expansion 

2.1 Basic PDE formulation 

Consider asset value processes X i: 1 < i < N, satisfying 



on a probability space {Q, J 7 , V} with filtration 6 [0,7"], T 6 IR + U oo. Here 

a : R N x [0, T] ->• M^' + is the volatility, /i : R N x [0, T] ->- R N is the drift, W Q is a standard 
Brownian motion under the risk-neutral measure Q and p : M. N x [0,T] — > R 7Vx7V is the 
correlation matrix, ie, 



A European option is characterised by its payout function G : M — > IR, which determines 
the amount G(Xt) its holder receives at time t = T. The arbitrage-free value of the option 
relative to the numeraire is then 



dXi(t) = /ii(X, t) dt + <Ti(X, t) dWf 



(2) 



(dWf,dWf)= Pij dt Vi,jel,...,N. 



(3) 




(4) 
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assuming: that standard technical conditions holcfl. Here G(-) is the absolute payoff at time 
T. By the Feynman-Kac theorem, V satisfies the parabolic PDE 



dV A dV A x d 2 V 

i=l ij=l 

on R N x [0, T] with boundary condition 

V(X,T)=g(X) VXgR^, (6) 

where for simplicity of notation we have used the relative payoff g(-) = G(-)/J\f(T). Equa- 
tion (jnj) naturally generalises to the Bermudan and Ratchet cases discussed later, which 
are modelled by the introduction of additional, intermediate conditions at a fixed, finite 
set of tenor times T%, . . . , T N . 

Assume now that p and a are constant and \x a function of t alone. Let £ be the 
covariance matrix, = UiPijUj for all 1 < i,j < N . Let Q G M ArxAr be the matrix of 
eigenvectors of £ and let the eigenvalues A, be sorted in descending order, ie, Ai > A2 > 
. . . > Ajv > 0. Then the coordinate transformation 



where 



leads to 



r = T-t , z = Q T \n{X) + f3(T), (7) 

AO-) = - E ^ (ir + l ^ {s) ds ) ' 1 - 1 - N > (8) 



N 



^-^E A ^ = ° V(r,,)G[0,T]xM-, (9) 
i=i * 

where V(t,X) = u(r,z) and 

u(z,0) = g(ex P [Qz]). (10) 

2.2 Taylor expansion 

Consider now u as a function also of the vector A = (Ai, . . . , Ajv) of eigenvalues. For any 
n > 1 and A G , we can define 5X = A — A and can formally write down the n-th order 
Taylor expansion at A as 

X du 



u(z,r;X) = u(z,T;X ) + J2 SX hQ^(z,r;X°) + ... 



N an 



h »n=l 

1 See, eg, [5]. 



+ ^ ^ • . . ■ • 5X in QXi U (z, r; A°) + O (\\SX\\ n+1 ) . (11) 
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The error term is justified for sufficient regularity of u. We can then choose suitable 
finite difference approximations A'' 1, "' , '" i '(fi) to d m u/d\i 1 . . -dXi m for all f < m < n and 
(ii, . . . , i m ) G {1, . . . , N} m . Hilber et al. [10] propose to use high order compact finite 
difference stencils introduced in [14"] . 

Choosing 5Xi as stepsize in direction i and denoting for each m by l m the lowest ap- 
proximation order of any we have 

( ' " (z, r; A ) = A^-^\u; z, r; A, A ) + O (\\6X\\ 1 ™) , (f 2) 



dX h ...d\i 

making explicit all arguments the finite difference approximation depends on. We addi- 
tionally set A°(w; z, r; A, A ) = u(z, r; A ). 

The finite difference approximation will contain the values u(z, r; A') for different values 
of A', which depend on A, A , and the finite difference formula itself. For all sensible finite 
difference approximations to derivatives of mixed order m, the number of non-zero elements 
of A' will bem<# plus the number of non-zeros of A . The computation of u(z, r; A') for 
a A' with k non-zero components can be accomplished by the solution of a fc-dimensional 
PDE of the form 

du 1 , d 2 u du I J^, , d 2 u 



i=l 1 i=l,A^0 2 

instead of the full A^-dimensional one. Insertion of ( Il2p in f lTT]) gives us 



A? 



u(z,r;A) = A°( M; z, r; A, A ) + . . . + ^ E ^ • . . . • 5A im A^>-^(«; r; A, A ) 

m=l ii,...,i m =l 

+ O (||5A|| min i<™<™('-+ m )) + (||5AH n+1 ) . (14) 

The overall approximation order is I = min{ra + l,mini< m < n (l m + m)} and the error is 
e = O (\\SX\\ l ) for smooth enough u. 

2.3 First-order, one-dimensional case 

A practically good choice of A and the number m of terms to include depends on the 
problem at hand. It is a common feature of processes with strong correlation that there 
is a dominant eigenvalue which is much larger than the rest of the spectrum. This is also 
the case for the model parameters illustrated in Fig. [T] in Section HI 

This motivates to expand up to first order, n — 1, around A = (Ai, 0, . . . , 0). Using a 
simple first-order forward finite difference approximation 

AW («; z, r; A, A ) = X° + 5X^) - u(z,r; X°) 

oXi 
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to the first derivative, where e; is the i-ih canonical basis vector, i > 1, we get a scheme 
with overall order / = min{2, 1 + 1} = 2. The corresponding error is of size 0(||<5A|| 2 ) = 
0{\\ + . . . + X 2 N ). To evaluate ( 114]) up to n — 1, we have to solve the one-dimensional PDE 



^-^1^-1 = (16) 



du 1 <9 2 w 
lfr~2 X ~d~4 

and the iV — 1 two-dimensional PDEs 

du 1 <9 2 m 1 d 2 u 

2 < i < N, and obtain the approximate solution 

N 

u^(z,t;X) = M(2,r;A ) + ^(u(z,r;A° + A i e i )-M(2,r;A )) (18) 

1=2 

AT 

= (2-JV)u(z,T;A )+5^u(z,r;A + A i e i ). (19) 

2=2 

Here, the superscript (1, 1) indicates that A has one non-zero element, and we are trun- 
cating the Taylor expansion after the first term. The largest dimension of any PDE to be 
solved is 1+1=2. 



3 Generalisations and relation to anchored ANOVA 

Anchored ANOVA-decompositions are used in [8] to obtain dimension- adaptive approxi- 
mations to option values expressed as integrals over high- dimensional spaces; [21 J point 
out a relation to the expansions from [20J by utilising the integral representation of the 
solution to the heat equation; I_QJ discusses these ideas jointly in the PDE context. 

We formulate the problem in slightly more general terms here as befits the applications 
later on. Consider the situation where Z = (Zi (£),..., Zjv(£))t>o is an N- dimensional 
Markov process and where the time t value u of a contingent claim is fully determined 
by Z(t) and the time to maturity T — t. We therefore write this value as u(Z(t), T — t). 
Define, for a set v = ■ ■ ■ ,i m } Q {l> - >^}> auxiliary processes Z v = (Zf,...,ZJ f ) 
which are "frozen" at the components with indices not in v, that is, Z?(t) = Zi(0) for all 
t > 0, i ^ v, and for all i e v we impose that the law of Z\ is identical to the law of Zi 
given that Zj(t) = Zj(0) for all j ^ v. Specifically, in the common case where Z is defined 
through a stochastic differential equation (SDE) of the type 

dZi(t) = i* z (Z(t),t) dt + a z {Z(t),t) dWt(t), (20) 

we define Z v by 

dZ v (t) = { VzA ZV (-t),t) dt + a ZA {Z v {t),t) dWiit) i e v, 
% I else, 
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which are constant in directions i G {1, . . . , N}\v, ie, Z"(t) = Zf(0) for those i. 

Now, define u v as the "solution of the problem with Z replaced by Z vn . We express this 
in such vague terms in order to include a wide class of applications. A particular example 
studied in [81 |2T], which is related to the valuation of European-style derivatives, is 

u{z, r) = u(z, T-t)= E[g(Z(T))\Z(t) = z], (22) 

where r = T — t is the time to maturity, g is the payoff function and the expectation is with 
respect to an explicitly known probability measure. Then, the "solution of the problem 
with Z replaced by Z vn can be written as 

u v (a; z\t) = E\g(Z%.)\Z?(t) = z { Vi G v; Z?(t) = a, V* £ v], (23) 

where, = Z^O), z v = (z^, . . . , Zi H ), ie, we anchor the solution at the initial value of this 
stochastic process. 

For optimal stopping problems, such as the Bermudan swaptions studied later, 

u(z,t) = su V E[g(Z(T))\Z(t) = z], (24) 
TeT 

u v (a;z v ,r) = sup E[g(Z v (T))\Z^(t) = z, Vi E v; ZV(t) = on Vi £ v], (25) 

where T, T v are suitable sets of stopping times. For path-dependent options, the process 
Z has to be set up to include the path-dependent quantity in order to bring it back into 
the assumed Markovian framework. In Section 15. 3[ we demonstrate this on the example 
of a ratchet floor. As the path-dependent quantity is reset at discrete time points, the 
corresponding component of Z is a jump-process instead of following an SDE of the form 

(EQD. 

The forward and backward PDEs for processes of the type (1201) are second order linear 
parabolic. To get from the PDE for (|22|) to the one for (|23|) . the coefficients of all derivatives 
in directions Xi are set to zero for i ^ v, as per (|2"T|) . 

The coordinates z in (171) were chosen specifically as the principal components of the 
covariance matrix of a diffusion process X, for constant \x and a, in which case fiz,i — 
and a\ i = \. A similar construction is used and analysed in [10]. The PDE satisfied by 
Uv in (E3D is (USD with A' set to 

X V = J2^- (26) 

As the solution u v , v = {ii, . . . ,i\ v \}, only depends on the sub-vector of coordinates 
z v = (z h ,. . . , Zi H ) non-trivially, 

u v (a;z v ,r) = u v (27) 

toCo 

= £(-l) H - H ^ (28) 

wCv 
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is well-defined and gives a suitable anchored ANOVA decomposition of u (see [5] and the 
references therein) as given in (JT]), where r is an additional argument of all terms. 

To re-iterate, the key is that Z'" changes only in \v\ dimensions and is constant in the 
remaining N — \v\ dimensions. This means that u v can be found from \v [-dimensional 
problems instead of the iV-dimensional one. 

The link between (JTJ) and f[T4"l) can now be established if we pick A = 0, set v = 
{ii,..., ik}, X v as in (|26|) . and, inductively (skipping z and r as argument of A for brevity), 

A^( W ;A,A°) = ^ T;A0 + ^ ) "^ r;A ° ) , (29) 
x,\°) = A (ifc+l) (A (il, ---'* fc) ; A, X v ) (30) 

(31) 



A(n,-«)( u; a, \ v + 5\ ik+1 e ik+1 ) - A^>-^(u; A, A") 



= ]T (-l)l«'-M-y 2(jT ;A«) l (32) 

U)CuU{fc+l} 

for ^ i m , 1 < ^ < m < k + 1, and otherwise. Only terms of mixed first order are 
present and absorb all higher order terms - see next paragraph. Then, the precise relation 
between ANOVA terms and the finite difference approximation to the Taylor expansion is 

u v {z;z v ,0) = u(z,0;\ v ), (33) 
u v {z-z v ,0) = A (il -"' ifc) (M;z,0;A,A°). (34) 

The relation between ANOVA and multi-variate Taylor expansions in the coordinates is 
discussed in [7]. The twist here is to apply the expansion in 5Xi instead of Zj. 

A further point to note is that if the above expansion is truncated to include terms up 
to \v\ — n < N, it is only of first order accurate in 5Aj. For relatively large 5Xi and smooth 
solutions, the inclusion of higher order Taylor terms in individual and mixed directions 
and higher order finite difference formulae like those proposed in [10] may be preferable. 
The extra cost is small as typically the dimensionality of PDEs involved will not increase. 
What distinguishes the above expansion from other finite difference approximations is that 
it is an exact decomposition, ie, if we include all terms up to degree N, we recover the 
exact solution irrespective of its smoothness. 

In a variation to ([T]), we can consider a decomposition, where in addition to the anchor 
a, all contributions may also depend on the first coordinate, 



TV N 



u(z) = Uq (a; 



zi) + u i (°i z r\ z i) + X] M li ( a ' Zl > Zi > z i) + 



i=2 i,i = 2 

i < j 

+ u ( 2 ) .. N (a;z 1 ;z 2j ... J z N ), (35) 
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and, generalising this from one to r > 1 coordinates, 

N-r 

u(z) = u { o\a;zi,...,z r ) + ^ ^ u { [^ _ >ip {a\ z u . . . , z r ; z h , . . . , z ip ). (36) 

p=i {h,...,i p } 

C{r + 1,...,N} 

Clearly, in relation to Section 12. 2\ this corresponds to using A = (Ai, 0, . . . , 0) and A = 
(Ai, A2, • • • , A r , 0, . . . , 0), resp., and adaptation of the finite difference formulae. 

The goal is to find a decomposition where the contributions decay fast with r + p, in 
order for the approximations 

s 

u {r ' s \z) = v%\a] z u . . . , Zr) + u C...,i p { a 'i z U---,z r ]Z il ,...,z ip ) (37) 

C {r + 1,...,N} 

for s < N—r to be accurate for small r+s. The approximation from Section l2T3"l corresponds 
to r = s = 1 . It is this case which will be used for the numerical studies in Section The 
effect of higher-dimensional terms in the cases r = 2, s — 1 and r = 1, s — 2 is illustrated 
in [21] for equity basket options, extending the case r = 1, s = 1 in [20J. The data there 
have in common with our set-up the presence of a dominant eigenvalue, such that the 
case (r, s) = (1, 2) gives a notable improvement over (1, 1) for arithmetic average basket 
options by capturing the second order terms in the small eigenvalues, while (r, s) = (2, 1) 
does not give a big accuracy gain. In our setting, we will see that additional parameter 
approximations are required, which can lead to errors comparable to the remainder terms 
in the expansion to order (1, 1). Consequently, the inclusion of higher order terms would 
not be worthwhile without a further extension of the method as sketched in Section 



4 Application to the LIBOR Market Model 

We now apply the PCA-ANOVA approach to practically relevant examples from interest 
rate markets. Forward rates will be assumed to follow the LIBOR Market Model (LMM), 
which is one of the most widely used models [31 El US] and the basis for a variety of 
extensions. The methods studied here have the potential to be applied to those as well. 
Our notation and definition of the LMM follows [5J. 

As traded product at time t consider a (zero-coupon) bond P(t,T), T > 0, that pays 1 
at time T > t. The forward LIBOR with fixing date T and payout date T' is then defined 
as the stochastic process L(-,T,T') : [0, T] x O — > R given by 

L{t,T,T) = wr - T . (38) 

For a fixed tenor structure 

= T x < . . . < T N < T N+l = T', (39) 
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the LMM now describes a finite number of forward rates 



Li(t) = L(t,T h T i+1 ) (40) 

for 1 < i < N, t e [0, Let Tj - = a for all 2 < j < N + 1. Practically important 
values for a are 0.25 and 0.5 for 3-month and 6-month LIBOR. The full dynamics for each 
Li(t), t <Ti, under the equivalent martingale measure Q n+1 , 1 < n < N, associated with 
choosing the bond P{T n+ i) as numeraire, are 



where 



for t < T, 




dL.it) = mWLiit) dt + ai{t)Li[t) dWf + ) (41) 

i = n (42) 



% > n 



N 



^ = E 1 T^^ ^^P^H^ (43) 

j'=max{fc:Tj.<t} 



1 + aLj(t) 



for t > T n+1 . 

Our model for the correlation structure is taken from [12J with 

Pij = exp(-#-j|) (44) 

for 1 < i,j < N and > and a constant volatility c = 0.2. The eigenvalues of the 
covariance matrix S decrease rapidly. Figure [1] demonstrates this for N — 11 and a typical 
set of parameters. The first eigenvalue Ai is roughly an order of magnitude larger then the 
second eigenvalue A2. This motivates the use of the first order, one-dimensional approach 
from Section 12.31 

We now choose the terminal bond P(T/v+i) as numeraire and combine the LIBOR 
dynamics in equations (|4"T|) - (|4"3"|) and our covariance structure with equation ([5]) to obtain 
a PDE satisfied by the value function of derivatives on the LIBOR curve. A complication 
arises in the transformation (J2J) to the heat equation (Q, as the drift term (3. in flH]) was 
assumed to depend only on r whereas with pi as in equation (|4"2"|) it also depends on Lj, 
j > i. 

To make the PCA approach directly applicable, we first approximate the drift term. A 
common approach in practice is to "freeze" the drift partially at its initial value by setting 

*<«> = - £ (45) 

j=i+l 3K ' 

This introduces an error in the drift that grows with Lj(t) — Lj(0), and the approximation 
can be expected to be reasonably accurate for moderate values of and a. We will 
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Figure 1: Eigenvalues A«, 1 < i < N, of £ for AT = 11, constant volatility c = 0.2, and 
= 0.02065 (+), = 0.0413 (*) and = 0.0826 (x), resp., see pi]) . 

confirm this numerically by comparing the PDE results to Monte Carlo estimates with and 
without drift approximation. A more accurate procedure is suggested in Section |6j 

A second point of consideration is that Lj(t) is only financially meaningful for t < Tj. 
In order not to have to change the underlying set of arguments of the value function, and 
hence the PCA, at every tenor time Tj, we consider "extended" LIBORs that are also 
defined for T, < t < T/y. In the case of constant p and a, a possible extension is obtained 
by demanding that Li(t) follows ( 141"]) for all < t < T/y. Note that the exact option value 
does not depend on Li(t) for t > Tj and is thus not affected by this extension. 

Applying the first order, one-dimensional PCA ANOVA approach from Section I2TB1 now 
leads to the approximate solution 

N 

u^ l \z,T]X) = u{z,t;\ ) + Y,Hz,t;\° + X i e i )-u{z,T;X )) (46) 

i=2 

N 

= (2-iV)- M 02,T;A o ) + 5>(z,r;A o + (47) 

i=2 

where A = (A 1 , 0, . . . , 0), 

r — T — t, z = Q T ln(L) + p( T ) (48) 

and 

N / 2 N T (ft) \ 

A(r) = -rX)Oii y-E u^L w* . 1 < z < iV- (49) 
i=i V fc=j+i fcl ; / 

Here, Q is the orthogonal matrix of eigenvectors of S, and A the vector of eigenvalues. 



11 



The initial condition for all PDEs is given by g(L) where g is the payoff at time T. The 
quantity of interest is u^(Z(0),T; A), where Z(0) = Q T ln(L(0)) + 0(T). 

5 Implementation and numerical results 

We chose two types of derivatives to test the flexibility and accuracy of the approach and 
benchmark against Monte Carlo results: 

• short- to long-running Bermudan swaptions, where the combination of high-dimensionality 
and early exercise presents challenges for PDE and MC methods; 

• a ratchet floor, where the path-dependency is conceptually straightforward to include 
in a MC solver and needs adaptation of the PDE solver. 

5.1 Implementation of PDE solvers 

To compute the approximate solution defined by equation f H6|) we need to numerically 
solve one- and two-dimensional PDEs of the type (I16p and (fT7|) . These are standard and 
we used the following approach. 

The computational domain is infinite in the ^-coordinates. To avoid the introduction 
of artificial boundary conditions necessary when localising the domain, for each coordinate 
Zi, we map the interval (— oo, oo) to (0, 1) via 

1 . , 1 

Ui = - arctan (7^ + q) + -, (50) 

7T 2 

with parameters 7, and q. Under a standard growth condition on the solution at infinity, 
the resulting PDE is fully specified without boundary conditions at G {0, 1}, because 
the resulting non-constant coefficients of the transformed diffusion-equation vanish suffi- 
ciently fast at the boundaries (see [20, 23J). For call type options such as the Bermudan 
swaption discussed below we apply a payout cutoff at a value g max = 1000, which does not 
significantly impact the computed option value. 

We consider an equidistant grid with J + 1 gridpoints along each axis, such that in 
original coordinates the mesh is denser in the interesting region, which depends on the 
LIBOR rates at r = 0. For instance, in the case Lj(0) = 0.1, which will be considered 
later, we choose 7$ and q such that LIBORs between 0.02 and 0.5 are mapped to the 
interval [0.1,0.9]. 

For the discretisation we use the Crank-Nicolson scheme with central spatial differ- 
ences. In the two-dimensional case, we combine this with an Alternating Direction Implicit 
(ADI) factorisation [T7], such that the resulting tridiagonal matrix systems can be solved 
efficiently in linear time (ie, proportional to the system size). As the coefficients of the 
PDEs are constant in time, an initial LU factorisation of the tridiagonal matrices gave 
significant further speed-up. 
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Depending on the derivative contract, there can be additional parameters and interface 
conditions to be taken into account. The two examples we considered are Bermudan 
swaptions, which offer early exercise rights at discrete points in time, and Ratchet floors, 
where a strike parameter is reset depending on LIBORs at tenor dates, which makes the 
payoff strongly path-dependent. We describe both in more detail in the next sections. 

5.2 Bermudan swaption 

A Bermudan (payer) swaption with strike price K can be exercised at one of a set of exercise 
dates {T ei , . . . , T ejv ,} C {T 1; . . . , T N }. Here, we consider an an example 3-month LIBOR, 
ie, a = 0.25 and T$ — a (i — 1) = 0.25 (i — 1) for 1 < % < N, and Bermudan swaptions which 
can be exercised yearly, ie, {T ei . . . , T 6jv , } = {Ti, T 5 , T 9 , . . . , Tat}, assuming that N — 1 can 
be divided by 4. If the Bermudan swaption is exercised at Tj then the holder receives a 
(payer) swaption with payout 

V Sw a P tion,i = max ■ ^[Lj(Tj) - K]/P(T h T j+1 ), O^j . (51) 

The value Vbs,x of a Bermudan swaption is thus determined by backward induction through 

Vbs,n = Vs W aption,N and 

V BS ,i = max(V BS ,i+i, V Swap tion,i) at i E {1, 5, . . . , N - 4}. (52) 

Between T, and T i+4 , % G {1,5,..., N — 4}, the value function Vbs,i satisfies the LMM PDE 
as discussed in Section m which we approximate by PCA and first order anchored ANOVA 
decomposition explained ibid. The interface condition fl52l) can easily be incorporated in 
the present PDE discretisation by evaluating fl52l) on the computational grid. 

As reference solutions for the PDE results we use Monte Carlo (MC) estimates. The 
numerical approximation of multi-dimensional American and Bermudan options by Monte 
Carlo methods is an area of active current research. We mention recent work on the com- 
putation of tight bounds via iteration approaches (eg [12]) and via pathwise optimisation 
(see [1]). Here, we use the well-established and popular primal-dual approach for exer- 
cise policy learning due to Andersen and Broadie (see |2J), which provides Monte Carlo 
estimates for a lower bound V^ c and upper bound V^ c + A to the true option value. 

In the simulations, we used Ni = 10 6 paths for learning the exercise policy (of type 
'exercise strategy 1' in the notation of (2]), N 2 = 10 7 paths to calculate the lower bound 
and N outer = 5000 and N inner = 1000 paths for the outer and inner MC runs to compute 
the upper bound. For the time discretization of the LMM SDEs we used the log-Euler 
scheme with Mmc — 5 time steps per interval of length a. In the tests with "frozen" 
drift (ie, lognormal LIBORs), the discretisation is exact for M > 1. For the PDE we used 
J = 601 grid points in every direction and the Crank-Nicolson scheme with Mpde = 10 
time steps per time interval of length a. The numerical parameters, summarised in Table 
(TJ were chosen such that the numerical error is small compared to the difference between 
PDE and MC solution and typically is of order 0.1% or less of the derivative value. 
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Parameter Value Description 



J 


601 


Number of grid points in each direction 


Mp DE 


10 


Timesteps per interval of length a in the PDE computation 


Nt 


10 6 


Number of MC paths for exercise policy learning 


N 2 


10 7 


Number of MC paths to calculate the lower bound 


outer 


5000 


Number of outer MC paths to calculate the upper bound 


N inner 


1000 


Number of inner MC paths to calculate the upper bound 


M MC 


5 


Timesteps per per interval of length a in the MC computation 



Table 1: Numerical parameters for the Bermudan swaption PDE and MC computation. 

Results for Bermudan swaptions at-the-money (ATM, K = 0.1) are shown in Table 
|2j The PDE results are compared to the values calculated by MC simulation with frozen 
and full drift, to disentangle the effects of the drift approximation on the one hand and 
the dimension reduction on the other. The model parameters chosen were <fi = 0.0413 and 
c = 0.2, with a flat initial LIBOR curve with L;(0) = 0.1, all identical to [12] to facilitate 
direct comparison. 



N 


Vmc 


a 


Vmc + A o 


°"Ao 


V MC + A o/2 


VpDE 






5 


1.75E-03 


9.02E-07 


1.75E-03 


0.00E-00 


1.75E-03 


1.76E-03 


1.18E-05 


0.68% 


11 


1.21E-02 


5.42E-06 


1.22E-02 


5.20E-06 


1.21E-02 


1.24E-02 


2.61E-04 


2.15% 


21 


3.05E-02 


1.14E-05 


3.15E-02 


4.65E-05 


3.10E-02 


3.14E-02 


4.03E-04 


1.30% 


41 


6.17E-02 


1.94E-05 


6.68E-02 


1.46E-04 


6.42E-02 


6.57E-02 


1.44E-03 


2.24% 


61 


8.23E-02 


2.29E-05 


9.10E-02 


2.16E-04 


8.67E-02 


9.04E-02 


3.77E-03 


4.35% 


81 


9.45E-02 


2.39E-05 


1.06E-01 


2.69E-04 


1.00E-01 


1.07E-01 


6.84E-03 


6.83% 


101 


1.01E-01 


2.38E-05 


1.14E-01 


3.82E-04 


1.08E-01 


1.18E-01 


1.02E-02 


9.45% 


5 


1.75E-03 


9.01E-07 


1.75E-03 


0.00E+00 


1.75E-03 


1.76E-03 


1.33E-05 


0.76% 


11 


1.21E-02 


5.43E-06 


1.22E-02 


5.16E-06 


1.21E-02 


1.24E-02 


2.57E-04 


2.11% 


21 


3.06E-02 


1.15E-05 


3.16E-02 


4.66E-05 


3.11E-02 


3.14E-02 


3.17E-04 


1.02% 


41 


6.19E-02 


1.98E-05 


6.77E-02 


1.63E-04 


6.48E-02 


6.57E-02 


8.62E-04 


1.33% 


61 


8.26E-02 


2.35E-05 


9.27E-02 


2.44E-04 


8.77E-02 


9.04E-02 


2.74E-03 


3.13% 


81 


9.49E-02 


2.45E-05 


1.07E-01 


2.85E-04 


1.01E-01 


1.07E-01 


6.13E-03 


6.07% 


101 


1.02E-01 


2.45E-05 


1.16E-01 


8.44E-04 


1.09E-01 


1.18E-01 


6.72E-02 


8.53% 



Table 2: PDE results for ATM (K = 0.10) Bermudan swaptions compared to frozen (top) 
and full (bottom) drift MC results. V^ [C and a are the lower MC bound and its estimated 
standard error. V^ c + A and <ta are the upper MC bound and the estimated standard 
error of the MC spread A . Vpde is the PDE result and the columns A a t, s and A re ; show 
the absolute and relative difference to the best MC estimate V^ IC + A /2. 

In these tests, the PDE method shows very good accuracy for up to N = 41: Vpde is 
above V^ c and below the upper MC bound in almost all cases for N = 21 or iV = 41. For 
lower N it is often slightly higher than the (in these cases fairly tight) upper bound, but 
the difference A abs to the middle value V^ c + A /2 - which is considered to be a better 
estimate for the true price than both the lower or upper bound in [12] - is always less than 
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2.6 basis points. For N = 61 — 101, the PDE values are still close to the MC values. They 
are above the upper MC bound for N — 81 and N = 101, though. Taking into account the 
relatively large difference between lower and upper M C bounds, this seems to indicate that 
the applicability of the first-order, one-dimensional version of the PCA-ANOVA approach 
reaches its limits (in this setting) for problems with N higher than 50 — 60 (as does the 
MC approach used here). 

In-the-money (ITM, K = 0.09) and out-of-the-money (OTM, K = 0.11) results are 
similar and are shown in Tables [TT] and [T2] in Appendix |A] Not surprisingly, the relative 
difference A re i = A a b s /(V^ IC + Aq/2) is typically smallest for ITM options and largest for 
OTM options. 

To assess the PCA-ANOVA approach under a range of market conditions, we present 
simulations for ATM Bermudan swaptions with differing parameters: Tables [3] and H] show 
the results for a smaller volatility of c = 0.1 and c = 0.15, resp., and Tables and [6] show 
the behaviour for weaker and stronger correlation between the LIBORs. As one would 
expect, the difference between PDE and MC solution decreases with lower volatility and 
increases with weaker correlation. Table [7] shows that the PCA approach also performs 
well for a lower initial LIBOR curve with Lj(0) = 0.02 for all 1 < i < N. 

Finally, Figures [2] and E] show how the value of the PDE approximation changes when 
we only consider the 1-dimensional PDE solution and the 2-dimensional PDE solutions 
associated with the k largest eigenvalues in in ( 146]) . Specifically, the case k = 1 is the 
one- dimensional approximation using the first principal component, and the case k = 2 
the standard two-dimensional PCA approximation. Evidently, it is in fact necessary to 
include the contributions from several of the largest eigenvalues to compute an accurate 
solution. At the same time, the solution levels out after including about 10 dimension. 
This is in line with the decay of the eigenvalues and the fact that the payoff in this case 
is almost parallel to the eigenvector of the first dimension. This suggests that for models 
where the eigenvalues of the covariance matrix decay fast enough - which includes many 
models in mathematical finance - the PCA-ANOVA PDE approach might be used with 
a fixed number k for a wide range of values N. This further reduces the computational 
effort without sacrificing significant accuracy. 

While considerable effort went into the efficient implementation of both the PDE and 
MC methods, there is still room for performance improvement, for instance on the algo- 
rithmic level and in the numerical parameter choices. Additionally, the computations were 
run on a number of different machines. Thus we only want to comment on approximate 
run times: For the PDE calculations, the computation times for N = 5, 11, 21 and 41 were 
on the order of 10, 60, 240 and 1200 seconds, resp., using Matlab on a AMD Phenom(tm) 
II X3 720 Processor (2.8 GHz) with 3.8 GB RAM. The run time is roughly quadratic in 
N, because the number of PDE solutions required to evaluate (l4"6"j) is N, and the expiry is 
T = N/4, so the number of (Crank-Nicolson) time steps with fixed step size, for a given 
PDE, is proportional to N. 

For the MC simulations the corresponding computation times were of order 300 + 
90, 1200 + 700, 4500 + 5000 and 18000 + 36000 seconds, resp. Here the first number is 
the computation time for the lower bound and the second number is the additional time 
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necessary to compute the upper bound. Again, the run time is quadratic in N, since N 
processes have to be simulated over N tenor dates. 

Despite the approximate nature of these computation times, it becomes clear that the 
PDE method is not only competitive time-wise, but indeed faster by a factor of 30 — 50 in 
our implementation. To get an optimal allocation of computational resources, one could 
also try to further optimise the relative size of the numerical parameters, such that the 
computation time is optimal for a given size of the combined discretisation error. By only 
halving the mesh size in the two directions of the computational grid one quadruples the 
computational time. There is no practical accuracy gain in bringing the discretisation error 
substantially below the error of the dimension reduction. A similar statement is true for 
the Monte Carlo estimators. A precise comparison of efficiency is therefore delicate. 

Both the PDE and MC methods used in this section are limited in their accuracy: the 
MC method uses a class of exercise strategies which generally does not include the optimal 
one; the PDE method employs a drift approximation and asymptotic expansion. For a 
wide range of N, the errors are comparable. Possible accuracy improvements on the basis 
of the decomposition are outlined in Sections |3] and |6j 



N 




a 


Vmc + A o 




V MC + A /2 


VpDE 






5 


8.78E-04 


4.29E-07 


8.78E-04 


0.00E+00 


8.78E-04 


8.82E-04 


3.32E-06 


0.38% 


11 


6.08E-03 


2.60E-06 


6.11E-03 


2.70E-06 


6.10E-03 


6.19E-03 


9.19E-05 


1.51% 


21 


1.53E-02 


5.57E-06 


1.59E-02 


2.45E-05 


1.56E-02 


1.56E-02 


2.87E-05 


0.18% 


41 


3.09E-02 


9.78E-06 


3.40E-02 


8.47E-05 


3.25E-02 


3.24E-02 


-1.40E-04 


-0.43% 


5 


8.78E-04 


4.29E-07 


8.78E-04 


0.00E+00 


8.78E-04 


8.82E-04 


3.32E-06 


0.38% 


11 


6.09E-03 


2.60E-06 


6.11E-03 


2.73E-06 


6.10E-03 


6.19E-03 


9.15E-05 


1.50% 


21 


1.53E-02 


5.54E-06 


1.58E-02 


2.18E-05 


1.56E-02 


1.56E-02 


5.50E-05 


0.35% 


41 


3.10E-02 


9.83E-06 


3.40E-02 


8.36E-05 


3.25E-02 


3.24E-02 


-1.44E-04 


-0.44% 



Table 3: PDE results for ATM (K = 0.1) Bermudan swaptions for a volatility of c = 0.1 
compared to frozen (top) and full (bottom) drift MC results. 



N 


*mc 


a 


Vmc + A o 


O-Aq 


V MC + A /2 


VpDE 






5 


1.31E-03 


6.59E-07 


1.31E-03 


0.00E+00 


1.31E-03 


1.32E-03 


8.05E-06 


0.61% 


11 


9.10E-03 


3.99E-06 


9.14E-03 


3.62E-06 


9.12E-03 


9.29E-03 


1.69E-04 


1.85% 


21 


2.30E-02 


8.48E-06 


2.37E-02 


3.41E-05 


2.33E-02 


2.35E-02 


1.74E-04 


0.75% 


41 


4.64E-02 


1.47E-05 


5.09E-02 


1.21E-04 


4.86E-02 


4.89E-02 


2.42E-04 


0.50% 


5 


1.31E-03 


6.59E-07 


1.31E-03 


0.00E+00 


1.31E-03 


1.32E-03 


9.04E-06 


0.69% 


11 


9.11E-03 


3.99E-06 


9.14E-03 


3.95E-06 


9.12E-03 


9.29E-03 


1.65E-04 


1.81% 


21 


2.30E-02 


8.50E-06 


2.38E-02 


3.63E-05 


2.34E-02 


2.35E-02 


1.12E-04 


0.48% 


41 


4.65E-02 


1.48E-05 


5.09E-02 


1.23E-04 


4.87E-02 


4.89E-02 


1.81E-04 


0.37% 



Table 4: PDE results for ATM (K = 0.1) Bermudan swaptions for a volatility of c = 0.15 
compared to frozen (top) and full (bottom) drift MC results. 
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N 


Vc 


a 


Vmc + A o 




Vmc + A o/2 


VpDE 






5 


1.75E-03 


9.02E-07 


1.75E-03 


0.00E+00 


1.75E-03 


1.76E-03 


1.54E-05 


0.88% 


11 


1.17E-02 


5.20E-06 


1.18E-02 


6.66E-06 


1.17E-02 


1.21E-02 


3.13E-04 


2.66% 


21 


2.85E-02 


1.05E-05 


2.95E-02 


4.39E-05 


2.90E-02 


2.96E-02 


6.59E-04 


2.27% 


41 


5.46E-02 


1.69E-05 


5.96E-02 


1.37E-04 


5.71E-02 


5.87E-02 


1.61E-03 


2.82% 


5 


1.75E-03 


9.02E-07 


1.75E-03 


0.00E+00 


1.75E-03 


1.76E-03 


1.56E-05 


0.89% 


11 


1.17E-02 


5.20E-06 


1.18E-02 


5.35E-06 


1.18E-02 


1.21E-02 


3.09E-04 


2.63% 


21 


2.85E-02 


1.05E-05 


2.96E-02 


4.74E-05 


2.91E-02 


2.96E-02 


5.52E-04 


1.90% 


41 


5.48E-02 


1.70E-05 


5.96E-02 


1.36E-04 


5.72E-02 


5.87E-02 


1.54E-03 


2.69% 



Table 5: PDE results for ATM (K = 0.1) Bermudan swaptions for weaker correlation with 
= 0.0413 x 2 = 0.0826 compared to frozen (top) and full (bottom) drift MC results. 



N 


Vmc 


a 


Vmc + A o 


°A 


V MC + A o/2 


VpDE 


A-abs 




5 


1.75E-03 


9.00E-07 


1.75E-03 


0.00E+00 


1.75E-03 


1.76E-03 


1.48E-05 


0.85% 


11 


1.23E-02 


5.54E-06 


1.23E-02 


5.35E-06 


1.23E-02 


1.26E-02 


2.43E-04 


1.98% 


21 


3.17E-02 


1.20E-05 


3.28E-02 


4.82E-05 


3.22E-02 


3.24E-02 


1.78E-04 


0.55% 


41 


6.64E-02 


2.11E-05 


7.17E-02 


1.57E-04 


6.90E-02 


6.99E-02 


8.26E-04 


1.20% 


5 


1.75E-03 


9.01E-07 


1.75E-03 


0.00E+00 


1.75E-03 


1.76E-03 


1.26E-05 


0.72% 


11 


1.23E-02 


5.56E-06 


1.24E-02 


5.84E-06 


1.24E-02 


1.26E-02 


2.14E-04 


1.74% 


21 


3.18E-02 


1.21E-05 


3.27E-02 


4.63E-05 


3.23E-02 


3.24E-02 


1.54E-04 


0.48% 


41 


6.66E-02 


2.16E-05 


7.28E-02 


1.70E-04 


6.97E-02 


6.99E-02 


1.59E-04 


0.23% 



Table 6: PDE results for ATM (K = 0.1) Bermudan swaptions for stronger correlation 
with = 0.0413/2 = 0.02065 compared to frozen (top) and full (bottom) drift MC results. 



N 


Vmc 


a 


V MC + A o 


°A 


V MC + A o/2 


VpDE 


A-abs 




5 


3.88E-04 


2.02E-07 


3.88E-04 


0.00E+00 


3.88E-04 


3.89E-04 


1.04E-06 


0.27% 


11 


2.86E-03 


1.31E-06 


2.87E-03 


1.37E-06 


2.87E-03 


2.90E-03 


3.78E-05 


1.32% 


21 


8.03E-03 


3.19E-06 


8.39E-03 


1.51E-05 


8.21E-03 


8.17E-03 


-3.72E-05 


-0.45% 


41 


2.00E-02 


7.34E-06 


2.32E-02 


7.56E-05 


2.16E-02 


2.08E-02 


-7.83E-04 


-3.63% 


5 


3.88E-04 


2.02E-07 


3.88E-04 


0.00E+00 


3.88E-04 


3.89E-04 


6.45E-07 


0.17% 


11 


2.86E-03 


1.31E-06 


2.88E-03 


1.46E-06 


2.87E-03 


2.90E-03 


3.74E-05 


1.30% 


21 


8.04E-03 


3.18E-06 


8.37E-03 


1.46E-05 


8.21E-03 


8.17E-03 


-3.31E-05 


-0.40% 


41 


2.01E-02 


7.36E-06 


2.30E-02 


7.41E-05 


2.15E-02 


2.08E-02 


-7.03E-04 


-3.27% 



Table 7: PDE results for ATM (K = 0.02) Bermudan swaptions in a flat initial LIBOR 
curve setting with Lj(0) = 0.02 VI < % < N compared to frozen (top) and full (bottom) 
drift MC results. 
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k 


VpDE 


A 


1 


5.37E-002 


-18.26% 


2 


6.08E-002 


-7.45% 


3 


6.29E-002 


-4.18% 


6 


6.47E-002 


-1.55% 


11 


6.52E-002 


-0.69% 


21 


6.55E-002 


-0.27% 


41 


6.57E-002 


0.00% 




Figure 2: Approximation Vpde to the value of the ATM Bermudan Swaption with N = 41 
when including only the terms up to i = k in (146)) . The column A shows the relative 
difference to the PDE result in the last row, which includes all N contributions. The solid 
horizontal lines (right) are the lower (primal) and upper (dual) MC bounds. 



•lO" 2 



k Vpde A 

I 2.85E-002 -9.40% 

2 3.04E-002 -3.34% q 

3 3.09E-002 -1.71% £ 
6 3.12E-002 -0.58% 

II 3.14E-002 -0.20% 
21 3.14E-002 0.00% 




5 10 15 

k 



Figure 3: Approximation to the value of the ATM Bermudan Swaption with N = 21 when 
including only the terms up to i = k in (j4U]) . The column A shows the relative difference 
to the PDE result in the last row, which includes all N contributions. The solid horizontal 
lines (right) are the lower (primal) and upper (dual) MC bounds. 
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5.3 Ratchet floor 



A Ratchet floor with strike price K\ and parameters a, b, c is a portfolio (sum) of floor- 
lets with payouts max{Ki — Lj(Tj),0} at T i+ i, where the strike prices Ki are recursively 
determined from the given initial strike K\ by 

Ki = max (aLi-i(Ti-i) + + c, 0) , i > 1, (53) 

see [16]. The (relative) price of the Ratchet floor for t = is given by 

N 

V RF (0) = J2 E iVFi,i( T i)/^( T i)\^o] , (54) 

i=l 

where Vp\^ is the value of the z-th floor let. Due to the linearity in the sum on the right- 
hand side of this equation it is sufficient to be able to calculate the price of a single floorlet. 
Without loss of generality we will thus focus on Vfi,n- 

The ratchet feature ( 153]) makes the problem path-dependent as the payoff depends on the 
values of all Li at different points in time. To solve the problem by a backward equation, 
we need to make it Markovian by including the strike dynamics with the evolution of the 
LIBORs, and to specify the value function as a function of all the above. As the strike 
changes discretely in time, the value function satisfies the standard LMM PDE between 
the Tj. At each tenor time Tj, the jump condition 

V F i ti (Ti-, K,Li,..., L N ) = V F i,i{Ti+, max (aLj_i + bK + c, 0), L x , . . . , L N ) 

holds. Here, Tj+ denotes the limit coming from larger t where the solution is already 
computed, and we use this to compute the solution just prior to Tj before the strike is 
updated. Details of the complete PDE model and its mathematical analysis can be found 
in [To|. 

We approximate the solution on a grid in ii'-direction, and compute the updated so- 
lution for each grid point via cubic spline interpolation for the corresponding value of 
Ki + \. This adds an extra dimension N + 1 to the problem, and effectively the one- 
dimensional ANOVA terms (corresponding to z\) now live on a two-dimensional grid, 
and the two-dimensional ANOVA terms (corresponding to Z\ and an additional z;) on a 
three-dimensional grid, formulaically (cf. Section [3]), 

N 

u^ l \K,z) = u^(K 1 ,Z(0);K,z 1 )+J2 u f\Ki,Z(0);K,z 1 ;z i ), (55) 

i=2 

where the superscript '2' stands for the variables K and Z\ and the '1' for the extra 
coordinate Z{ in the expansion. It is conceivable to use an expansion in direction K with 
anchor Ki of the form 

N 

u^\K, z) = u^iK,, Z(0); z x ) + £ uf\K 1} Z(0); z x] + v$ +l {K u Z(0); z x \ K), (56) 

i=2 
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Parameter Value Description 



J 
K 

M PDE 

N x 

M MC 



401 
81 
10 

10 7 
10 



Number of grid points in each LIBOR direction 
Number of grid points in strike price direction 
Timesteps per interval of length a in the PDE computation 
Number of MC paths 

Timesteps per per interval of length a in the MC computation 



Table 8: Numerical parameters for the Ratchet Floor PDE and MC computation. 

where the first superscript '1' stands for z\ and the second '1' for the expansion in one 
coordinate z% or K. This could work well especially if a, b > such that K and have 
positive correlation. We do not pursue this here, however. Due to the smoothness of the 
solution in the indirection and the higher order of the spline interpolation, a relatively 
coarse mesh is sufficient. Specifically, we use 21, 41 and 21 spline nodes for the intervals 
[0,0.05], [0.05,0.15] and [0.15, K max ), resp., with K max = 0.5, to give a total of 81 nodes. 
In comparison, we use J = 401 grid points in the Zj-directions and the Crank-Nicolson 
scheme with M — 10 time steps per interval of length a = 0.25 between tenors. The 
numerical parameters are again chosen such that the numerical error is small compared to 
the difference between PDE and MC solution and typically is of order 0.1% or less of the 
derivative value. The model parameters are identical to those for the Bermudan swaption 
in Section l5~2l 

We again use MC estimates as reference solutions. For a path-dependent option without 
early-exercise features like the Ratchet floor we can use a straightforward MC calculation. 
In the tests, we sampled Ni = 10 7 paths; for the time discretization, we used the log-Euler 
scheme with M = 10 time steps per interval of length a. 

The numerical results for the Ratchet floors are shown in Tables [9] and [10] for N = 
5, 11, 21. We consider three different configurations (a, b, c) and three different strike prices 
K = 0.1,0.11,0.09 (ATM, ITM and OTM, resp.). The absolute difference between the 
PDE and MC solution is never more than 1.14 basis points and the relative difference is 
below 1% in all but one case. 

We again report approximate computation times with the same caution as in the previous 
section. For N = 5, 11 and 21 the PDE run times were of order 800, 5000 and 20000 
seconds, resp., on a AMD Phenom(tm) II X3 720 Processor (2.8 GHz) with 3.8 GB RAM, 
while the MC run times were of order 600, 3000 and 11000 seconds, resp. The addition of 
an extra dimension for K in the PDE computation resulted in the MC computation to be 
faster by a factor of 1.3 to 2. The MC simulation also permits the computation of values 
for multiple parameters (a, b, c) in parallel, with only a small increase in computation time. 
Given the fast decay of the correction terms in the ANOVA decomposition, as illustrated in 
Figures [2] and [3] for the Bermudan swaption, it would be possible to compute only the first 
k ^ N ANOVA terms without significant loss of accuracy, which brings the computational 
times for the PDE in the range of, or below, the MC ones. 
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Table 9: PDE results for Ratchet floors compared to frozen drift MC results. The column 
Vpde shows the computed PDE value. Columns A a f, s and A re i show the absolute and 
relative difference to the MC estimate Vmc, resp. 
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Table 10: PDE results for Ratchet floors compared to full drift MC results. The column 
Vpde shows the computed PDE value. Columns A a f, s and A re i show the absolute and 
relative difference to the MC estimate Vmc, resp. 
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6 Conclusion and outlook 



The results presented in this article demonstrate the practical applicability of a systematic 
expansion approach to the LIBOR Market Model. We were able to compute values for 
Bermudan swaptions and Ratchet floors which showed a very good match to the Monte 
Carlo benchmark for up to N = 40 — 60 quarterly LIBORs under a range of market condi- 
tions. The deterioration in the accuracy for larger N is less an effect of the dimensionality 
as such as of the increased time to maturity. Based on the structure of both the problem 
and the method, and also encouraged by results reported in [21] for vanilla basket options, 
we are optimistic that the inclusion of higher order terms would give a practically sufficient 
fit also for large values of N if those were needed. The run times were of the same order of 
magnitude as the MC run times for the path-dependent example and substantially smaller 
for the early exercise case. 

Further work is needed for the analysis. Section[2]motivates the PCA-ANOVA approach 
via Taylor expansion, but it does not fully prove its applicability or show in which cases 
it cannot be used. In particular, the size of the coefficients in the Taylor and ANOVA 
expansions depends on the smoothness of the solution and it is ongoing work to derive 
error bounds. Encouragingly, a closer look at the raw data going into Figures [2] and [3] 
reveals that the correction terms are indeed strictly decreasing. This could be used as the 
(heuristic) basis for dimension adaptivity. 

A practically important extension is to treat variable coefficients in the PDE accurately 
and systematically. In the model studied here, the covariance matrix was assumed constant 
and the non-constant drift was approximated by a constant one. A more general approach 
to variable coefficients would be to "freeze" only a subset of the covariance and drift 
components as is required for the anchored ANOVA. We expect this to give higher order 
accuracy in T. This will allow us to use more complex volatility and correlation structures 
than the ones described in Section 0J 

Overall, we believe that the approach discussed here can be developed into an extremely 
powerful and versatile framework for the approximation of high-dimensional problems. It 
is not inherently restricted to the LIBOR market or Mathematical Finance problems in 
general and we expect it to perform well across a wide range of problems with suitable 
correlation structures. 
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A Further results 
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Table 11: PDE results for ITM (K = 0.09) Bermudan swaptions compared to frozen (top) 
and full (bottom) drift MC results. 
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Table 12: PDE results for OTM (K = 0.11) Bermudan swaptions compared to frozen (top) 
and full (bottom) drift MC results. 
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This figure "lambdall.png" is available in "png" format from: 



http://arxiv.org/ps/1209.1909vl 



